home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_6 / a6_1.m next >
Encoding:
Text File  |  1994-06-05  |  2.4 KB  |  106 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 6.1 (Differentiation Using Limits).
  9. % Section    6.1,    Approximating the Derivative, Page 326
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program finds a numerical approximation for f'(x0).
  13.  
  14. % The method employed is the limit process.
  15.  
  16. %    The function f(x) is stored in the M-file  f.m.
  17. % function z = f(x)
  18. % z = cos(x);
  19.  
  20. delete f.m
  21. diary f.m; disp('function z = f(x)');...
  22.            disp('z = cos(x);');...
  23. diary off;
  24.  
  25. f(0); % Remark. f.m and difflim.m are used for Algorithm 6.1
  26.  
  27. pause % Press any key to continue.
  28.  
  29. clc;
  30. %    Place abscissa for differentiation in x0.
  31.  
  32. % Place the endpoints interval [a,b] about x0 in a and b.
  33.  
  34. x0 = 0.8;
  35.  
  36. a  = 0;
  37.  
  38. b  = pi/2;
  39.  
  40. y0 = f(x0);
  41.  
  42. pause    % Press any key to graph the function.
  43.  
  44. clc; clg;
  45. h  = (b-a)/100;
  46. Xs = a:h:b;
  47. Ys = f(Xs);
  48. a = 0;
  49. b = 1.6;
  50. c = 0;
  51. d = 1.1;
  52. axis([a b c d]);...
  53. plot(x0,y0,'or',Xs,Ys,'g');...
  54. hold on;...
  55. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  56. xlabel('x');...
  57. ylabel('y');...
  58. title('y = f(x) and the given point');...
  59. grid;...
  60. axis;...
  61. hold off;...
  62. shg; pause    % Press any key to continue.
  63.  
  64. clc;
  65. %    Place abscissa for differentiation in x0.
  66.  
  67. % Place the tolerance in  toler.
  68.  
  69. x0 = 0.8;
  70.  
  71. toler = 1e-12;
  72.  
  73. [H,D,E,n] = difflim('f',x0,toler);
  74.  
  75. pause    % Press any key to see the approximate derivative.
  76.  
  77. X1 = [a b];
  78. C = [D(n) f(x0)];
  79. Y1 = polyval(C,X1-x0);
  80. axis([a b c d]);...
  81. plot(x0,y0,'or',Xs,Ys,'-g',X1,Y1,'-r');...
  82. hold on;...
  83. plot([-0.05 1.6],[0 0],'b',[0 0],[-0.05 1.05],'b');...
  84. xlabel('x');...
  85. ylabel('y');...
  86. title('The tangent line has slope m = f`(x0).');...
  87. grid;...
  88. axis;...
  89. hold off;...
  90. shg; pause    % Press any key to continue.
  91.  
  92. points = [H;D];
  93. Mx0 = 'Approximating the derivative by finding a limit.';
  94. Mx1 = 'Numerical approximations for ';
  95. Mx2 =  ['f`(',num2str(x0),')'];
  96. Mx3 = [Mx1,Mx2];
  97. Mx4 = 'The approximate value of ';
  98. Mx5 = [Mx4,Mx2,' = '];
  99. Mx6 = 'An estimate for the error bound is:';
  100. Mx7 = ' approx. derivative   +- error bound';
  101. clc,echo off,diary output,...
  102. disp(''),disp(Mx0),disp(''),disp(Mx3),...
  103. disp('     h                  D(k)'),disp(points'),...
  104. disp(''),disp(Mx5),disp(D(n)),...
  105. disp(Mx6),disp(Mx7),disp([D(n) E(n)]),diary off,echo on
  106.